double setDiameters(int npart, double sigma, double *d)
{
	double vol=0.0;
	/*Calculamos el diametro de las particulas*/
	/*(distribucion gaussiana media 1 variancia sigma)*/
	for(int i=0;i<npart;i++)
	{
		double a = RandomGaussianNumber();
		d[i]=sigma*RandomGaussianNumber()+1.0;
		vol+=POW3(d[i]);
	}
	return vol;
}

void setPositions(int npart, double **pos, double *d, double rmax, double ho)
{
	double u1, u2, r;
	bool ok;
	/*Calculamos las posiciones*/
	r=0.0;
	for(int i=0;i<npart;i++){
		ok=false;
		while(!ok){
			ok=true;
			u1=RandomNumber();
			u2=RandomNumber();
			pos[i][0]=rmax*u1*cos(2.*pi*u2);
			pos[i][1]=rmax*u1*sin(2.*pi*u2);
			u1=RandomNumber();
			pos[i][2]=(ho-2.0*d[i])*u1+d[i];
			for(int j=0;j<i;j++){
				//if(i==286 && j==285)
				//{
				//	double d1=(pos[i][0]-pos[j][0]);
				//	double d2=(pos[i][1]-pos[j][1]);
				//	double d3=(pos[i][2]-pos[j][2]);
				//	double re=sqrt(d1*d1+d2*d2+d3*d3);
				//	cout<<"a"<<endl;
				//}
				r=sqrt((pos[i][0]-pos[j][0])*(pos[i][0]-pos[j][0])+
					(pos[i][1]-pos[j][1])*(pos[i][1]-pos[j][1])+
					(pos[i][2]-pos[j][2])*(pos[i][2]-pos[j][2]));
				if (r<1.1*(d[i]+d[j])/2.){
					ok=false;
				}
			}
		}
	}
}

void calc_force (int npart, double **pos, double *d, double **ft, double wt, double vwd)
{
	/*componentes de fuerza repulsiva de corto alcance*/
	double frT,fr[3];
	/*componentes de posiciones relativas interparticula*/
	double r,r2, rm4,dr[3];
	double fm[3],dij3;
	/*Coseno(Theta) y su cuadrado*/
	double Co, Co2;
	/*Variables auxiliares*/
	double ang1, ang2, t1;

	/*Inicializamos la fuerza total*/
	for(int i=0;i<npart;i++){
		for(int k=0;k<3;k++){
			ft[i][k]=0.0;
		}
	}

	for(int i=0;i<npart-1;i++){
		for(int j=i+1;j<npart;j++){
			dij3=POW3(d[i])*POW3(d[j]);
			r2=0.0;
			/*Calculamos distancias*/
			for(int k=0;k<3;k++){
				dr[k]=pos[j][k]-pos[i][k];
				r2+=dr[k]*dr[k];
				fr[k]=0.;
			}
			r=sqrt(r2);
			rm4=1./(r2*r2);
			double rrr=sqrt((pos[i][0]-pos[j][0])*(pos[i][0]-pos[j][0])+
					(pos[i][1]-pos[j][1])*(pos[i][1]-pos[j][1])+
					(pos[i][2]-pos[j][2])*(pos[i][2]-pos[j][2]));

			/*Calculamos las fuerzas repulsivas*/
			if (r<(d[i]+d[j])){
				frT=-exp(-100*(2*r/(d[i]+d[j])-1.));
				for(int k=0;k<3;k++){
					fr[k]=frT*dr[k]/r;
				}			
			}
			/*Calculamos las fuerzas magneticas*/
			Co=dr[2]/r;
			Co2=Co*Co;
			t1=dij3*rm4;
			ang1=t1*(5*Co2-1.0);
			ang2=t1*(5*Co2-3.0);
			fm[0] = ang1*(dr[0]/r);
			fm[1] = ang1*(dr[1]/r);
			fm[2] = ang2*(dr[2]/r);
			/*Calculamos las fuerzas totales*/
			for(int k=0;k<3;k++){
				t1=fm[k]+fr[k];
				ft[i][k]+=t1;
				ft[j][k]-=t1;
			}
		}

		/*Calculamos las componentes de las fuerzas hidrodinamicas dependientes de la velocidad laminar u0*/
		/*(Por el funcionamiento del bucle falta la ultima particula)*/
		ft[i][0]+=3.0*pos[i][2]*vwd*d[i]/wt/wt*(1.-(pos[i][2]/wt))*pos[i][0];
		ft[i][1]+=3.0*pos[i][2]*vwd*d[i]/wt/wt*(1.-(pos[i][2]/wt))*pos[i][1];
		ft[i][2]-=3.0*pos[i][2]*pos[i][2]*vwd*d[i]/wt/wt*(1.0-2.0/3.0*(pos[i][2]/wt));


		/*Fuerzas de las particulas con las paredes*/
		/*(Por el funcionamiento del bucle falta la ultima particula)*/
		if (pos[i][2]<=d[i]){
			ft[i][2]+=exp(-100.0*(pos[i][2]/d[i]-0.5));
		}
		if (pos[i][2]>= (wt-d[i])){
			ft[i][2]-=exp(-100.0*((wt-pos[i][2])/d[i]-0.5));
		}
	}

	/*Incluimos la fuerza hidrodinamica y con las paredes de la particula que faltaba*/
	ft[npart-1][0]+=3.0*pos[npart-1][2]*vwd*d[npart-1]/wt/wt*(1.0-(pos[npart-1][2]/wt))*pos[npart-1][0];
	ft[npart-1][1]+=3.0*pos[npart-1][2]*vwd*d[npart-1]/wt/wt*(1.0-(pos[npart-1][2]/wt))*pos[npart-1][1];
	ft[npart-1][2]-=3.0*pos[npart-1][2]*pos[npart-1][2]*vwd*d[npart-1]/wt/wt*(1.-2./3.*(pos[npart-1][2]/wt));

	if (pos[npart-1][2]<=d[npart-1]){
		ft[npart-1][2]+=exp(-100.0*(pos[npart-1][2]/d[npart-1]-0.5));
	}

	if (pos[npart-1][2]>=(wt-d[npart-1])){
		ft[npart-1][2]-=exp(-100.0*((wt-pos[npart-1][2])/d[npart-1]-0.5));
	}

}

void move_balls(int npart, double **pos, double **ft, double *d, double deltat,double wt)
{
	/*Consideramos particulas inerciales --> D2r/Dt2=0 y por tanto las sumas de las fuerzas es igual a 0*/
	/*De la componente de la fuerza viscosa despejamos dr/dt e integramos para obtener las posiciones*/
	double deltax,deltay;

	for(int i=0;i<npart;i++){
		if (pos[i][2]<=1.1*d[i]/2.0 || pos[i][2]>=(wt-(1.1*d[i]/2.0))){
			/*Si las particulas estan en contacto con las paredes no es posible el desplazamiento
			de la misma en el plano xy --> deltax=deltay=0*/
			deltax=0.0;
			deltay=0.0;
		}else{
			deltax=deltat*ft[i][0]/d[i];
			deltay=deltat*ft[i][1]/d[i];
		}
		pos[i][0]+=deltax;
		pos[i][1]+=deltay;
		pos[i][2]+=deltat*ft[i][2]/d[i];
	}
}

void calc_energy(int npart, double **pos, double *d, double ho, double wd,double t)
{
	double def,nn;
	double ut, a;
	double r, r2, dr[3];
	a=sqrt(6.0);
	ut=0.0;
	nn=0.;
	def=(ho-wd)/ho;
	for(int i=0;i<npart-1;i++){
		for(int j=i+1;j<npart;j++){
			r2=0.0;
			for(int k=0;k<3;k++){
				dr[k]=pos[j][k]-pos[i][k];
				r2+=dr[k]*dr[k];
			}
			r=sqrt(r2);

			if (r<1.25*(d[i]+d[j])/2.0){
				nn++;
			}
			ut+=POW3(d[i]*d[j])*(1./3.-(dr[2]/r)*(dr[2]/r))/POW3(r);
		}
	}
	ut/=double(npart);
	nn*=2.0/double(npart);
	printf("Energia: %5.5f %5.5f %5.5f %5.5f %5.5f\n",t,wd,def,ut,nn);
}

double distancia2(double **pos, int i, int j)
{
	double r2, dr;
	r2=0.0;
	for(int k=0;k<3;k++){
		dr=pos[j][k]-pos[i][k];
		r2+=dr*dr;
	}
	return dr;
}